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Abstract 

A recent calculation of the nuclear energy density functional from chiral two- and 
three-nucleon forces is extended to the isovector terms pertaining to different proton and 
neutron densities. An improved density- matrix expansion is adapted to the situation of 
small isospin-asymmetries and used to calculate in the Hartree-Fock approximation the 

^S| ' density-dependent strength functions associated with the isovector terms. The two-body 

interaction comprises of long-range multi-pion exchange contributions and a set of contact 
terms contributing up to fourth power in momenta. In addition, the leading order chiral 
three-nucleon interaction is employed with its parameters fixed in computations of nuclear 
C^ ■ few-body systems. With this input one finds for the asymmetry energy of nuclear matter 

the value A{pq) ~ 26.5 MeV, compatible with existing semi-empirical determinations. The 
strength functions of the isovector surface and spin-orbit coupling terms come out much 

CN ' smaller than those of the analogous isoscalar coupling terms and in the relevant density 

range one finds agreement with phenomenological Skyrme forces. The specific isospin- 
and density-dependences arising from the chiral two- and three-nucleon interactions can 
be explored and tested in neutron-rich systems. 
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> : 1 Introduction 

z^ • The nuclear energy density functional approach is the many-body method of choice in order 
\0 . to calculate the properties of medium- mass and heavy nuclei in a systematic manner [H [2]. 
ff^ [ Parameterized non-relativistic Skyrme functionals [31 H] as well as relativistic mean-field models 
^ I O E] have been widely and successfully used for such nuclear structure calculations. In a 
complementary approach one attempts to constrain the analytical form of the functional and the 
values of its couplings from many-body perturbation theory and the underlying two- and three- 
nucleon interaction. Switching from conventional hard-core NN-potentials to low-momentum 
interactions [71 [8] is essential in this respect, because the nuclear many-body problem formulated 
in terms of the latter becomes significantly more perturbative. 

In many-body perturbation theory the contributions to the energy are written in terms of 
density-matrices convoluted with the finite-range interaction kernels, and are therefore highly 
non-local in both space and time. In order to make such functionals numerically tractable in 
heavy open-shell nuclei it is necessary to develop simplified approximations for these functionals 
in terms of local densities and currents. In such a construction the density-matrix expansion 
comes prominently into play as it removes the non-local character of the exchange (Fock) contri- 
bution to the energy by mapping it onto a generalized Skyrme functional with density-dependent 
couplings. For some time the prototype for that has been the density-matrix expansion of Negele 
and Vautherin [9], but recently Gebremariam, Duguet and Bogner [10] have developed an im- 
proved version for spin-unsaturated nuclei. They have demonstrated that phase-space averaging 
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techniques allow for a consistent expansion of both the spin-independent (scalar) part as well 
as the spin-dependent (vector) part of the density-matrix. 

By applying these new techniques a microscopically constrained nuclear energy density func- 
tional has been derived from the chiral NN-potential at next-to-next-to-leading order (N^LO) in 
ref. [11] by Gebremariam, Bogner and Duguet. These authors have proposed that the density- 
dependent couplings associated with the pion-exchange interactions should be added to a stan- 
dard Skyrme functional (with several adjustable parameters). In the sequel it has been demon- 
strated in ref. [12] that this new energy density functional gives numerically stable results and 
that it exhibits a small but systematic reduction of the x^'^^G^isure compared to standard 
Skyrme functionals (without any pion-exchange terms). 

In the recent work [13] the calculation of the nuclear energy density functional has been 
continued and extended with improved (chiral) two- and three- nucleon interactions as input. 
For the two-body interaction the N^LO chiral NN-potential has been used in ref. p3] . It consists 
of long-range multi-pion exchange terms and two dozen low-energy constants which param- 
eterize the short- distance part of the NN- interaction. The actual calculation in ref.[T3] has 
been performed with the version N'^LOW developed in ref. [13] by lowering the cut-off scale to 
A = 414 MeV. This value coincides with the resolution scale below which evolved low-momentum 
NN-potentials become nearly model-independent and exhibit desirable convergence properties 
in perturbative many-body calculations [TJEllin]. The (low-momentum) two-body interaction 
N^LOW has been supplemented in ref. [13] by the leading order (N^LO) chiral three-nucleon in- 
teraction with its parameters ce, cd and 01,3,4 determined in computations of nuclear few-body 
systems [151 [16]. With this input the nuclear energy density functional has been derived to first 
order in many-body perturbation theory, i.e. in the Hartree-Fock approximation. For the effec- 
tive nucleon mass M*{p) and the strength functions -Fv(p) and Fso{p) of the (isoscalar) surface 
and spin-orbit coupling terms reasonable agreement with results of phenomenological Skyrme 
forces has been found (in the relevant density range). However, as indicated in particular by the 
nuclear matter equation of state E{p), an improved description of the energy density functional 
requires at least the treatment of the two-nucleon interaction to second order in many-body 
perturbation theory. 

The purpose of the present paper is to extend the calculation of the nuclear energy density 
functional in ref. [13] to isospin-asymmetric many-nucleon systems with different proton and 
neutron densities. The additional isovector terms play an important role in the description 
of long chains of stable isotopes and for nuclei far from stability. Our paper is organized as 
follows. In section 2 we recall the improved density-matrix expansion of Gebremariam, Duguet 
and Bogner [TD] whose Fourier transform to momentum space provides the adequate technical 
tool to calculate the nuclear energy density functional in a diagrammatic framework. In section 
3 we present the two-body contributions to the various density-dependent strength functions 
A{p), Gr{p), Gd{p), Gso{p) and Gj{p), separately for the finite-range pion-exchange and the 
zero-range contact interactions. Section 4 comprises the corresponding analytical expressions 
for the three-body contributions grouped into contact (ce), Ivr-exchange (c^) and 27r-exchange 
(01,3,4) terms. Finally, we discuss in section 5 our numerical results and add some concluding 
remarks. 

2 Density-matrix expansion and isovector part of energy 
density functional 

The starting point for the construction of an explicit nuclear energy density functional is the 
bilocal density-matrix as given by a sum over the orbitals occupied by protons and neutrons: 



X]a^p"2(^— 0'/2)^I/p"2^(^+ 0/2). According to Gebremariam, Duguet and Bogner [TO] it can 
be expanded in relative and center-of-mass coordinates, a and r, with expansion coefficients 
determined by local proton and neutron densities. These are the particle densities Pp,n{'f^), the 
kinetic energy densities Tp^ni''^) and the spin-orbit densities Jp^nif) (for definitions in terms of the 
orbitals ^^°'}^{r)^ see section 2 in ref.[I3j). The Fourier transform of the expanded density- matrix 
with respect to both coordinates defines in momentum space a medium insertion: 
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for the inhomogeneous isospin-asymmetric many-nucleon system. Here, T3 denotes the third 
Pauli isospin-matrix and we have displayed only the (relevant) terms proportional to differences 
of proton and neutron densities: Pp — pn, Tp — Tn, Jp — Jn- The local Fermi momenta kp^njif) 
are related to the (particle) densities in the usual way: pp = /Cp/Svr^, p„ = /c^/Stt^, p = Pp + Pn = 
2/cj/37r^. When working to quadratic order in deviations from isospin symmetry (i.e. proton- 
neutron differences) it is sufficient to use an average Fermi momentum kj in the prefactors of 
Tp T„ anci Jp jfi- 

Up to second order in proton-neutron differences and spatial gradients the isovector part of 
the nuclear energy density functional takes the form: 

^iv[Pp, Pn, Tp, r„, Jp, Jn] = -{Pp- Pnf ^(p) + -{.Tp - rn){Pp - Pn) Gt{p) 

+ (Vpp - Vpnf Gsj{p) + (Vpp - Vpn) ■ {Jp - Jn) G so{p) + {Jp " Jn)' Gj{p) . (2) 

Here, A{p) is the interacting part of the asymmetry energy of (homogeneous) nuclear matter. 
The non-interacting (kinetic energy) contribution v4kin(p) = k'j/QM to the asymmetry energy 
is included in the nuclear energy density functional through the kinetic energy density term, 
•^kin = (Tp + Tn)/2M, with M = 939 MeV the (free) nucleon mass. The strength function G\7{p) 
of the isovector surface term (Vpp — Vpn)"^ has the decomposition: 

G^{p) = j-GAP) + G,{p), (3) 

where Gd{p) comprises all those contributions for which the (Vpp — Vp^)^ factor originates 
directly from the momentum dependence of the interactions in an expansion up to order q'^. 
The Fourier transformation in eq.(l) converts this factor q"^ into (Vkp — Vkn)"^ ~ (Vpp — 
Vpn)^(7r/A;/)^. The second last term {'Vpp—'Vpn)-{Jp—Jn) Gso{p) in eq.(2) describes the isovector 
spin-orbit interaction in nuclei. Depending on the sign and size of its strength function Gso{p) 
the spin-orbit potentials for protons and neutrons are differently composed from the gradients 
of the local proton and neutron densities. 

3 Two-body contributions 

In this section the two-body contributions to the various strength functions A{p), Gr{p), Gd{p), 
Gso{p) and Gj{p) are worked out. We follow closely section 3 in ref. [13] where the input two-body 
interaction, the chiral nucleon- nucleon potential N^LOW [14], has been described in sufficient 



detail. In the (first-order) Hartree-Fock approximation the finite-range multi-pion exchange 
interactions lead in combination with the density-matrix expansion (i.e. by employing the 
product of two medium insertions T{pi, q) T{p2, —q)) to the following two-body contributions 
to the strength functions: 
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with the (isoscalar minus isovector) combination of the central, spin-spin and tensor NN- 
potentials in momentum space: 

U{q) = Vciq) - Wciq) + 5Vs{q) - 3Ws{q) + q'Vriq) - q'WM . 
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In the integrands of eqs.(4,9) the momentum transfer variable q is to be set to g = 2xkf. The 
double-prime in eq.(7) denotes a second derivative and we have given the numerical value for 
Gd{p) resulting from the (negative) curvature of the isovector central potential Wc{q) shown in 
Fig. 1 of ref.|l3j. 

In addition there are the two-body contributions from the zero-range contact potential of 
the chiral NN-interaction N^LOW. The corresponding expression in momentum space includes 
constant, quadratic, and quartic terms in momenta and it can be found in section 2.2 of ref. |17]. 
The Hartree-Fock contributions from the NN-contact potential to the strength functions read: 
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The 24 low-energy constants Cs,t, Gj and Dj are determined (at the cut-off scale of A = 
414 MeV) in fits to empirical NN-phase shifts and deuteron properties \l^. Their numerical 
values have been extracted from the pertinent NN-scattering code and are listed in section 3 
of ref. [13]. Let us mention that the contributions proportional to Cs,t and Cj in eqs. (10-14) 
have also been worked out in appendix B of ref. [TT] and we find agreement with their results. 
The terms proportional to Dj as well as the master formulas eqs. (4-9) for the finite-range 
contributions are new. 





Figure 1: Three-body diagrams related to the contact (ce) and Ivr-exchange (cd) component 
of the chiral three-nucleon interaction. The short double-hne symbohzes the medium insertion 
r{p,q) for inhomogeneous isospin-asymmetric nuclear matter. 

4 Three-body contributions 

In this section the three-body contributions to the strength functions A{p), Gr{p), Gd{p), Gso{p) 
and Gj{p) are worked out. We employ the leading order chiral three-nucleon interaction [TB] 
which consists of a contact piece (with parameter ce), a Ivr-exchange component (with parameter 
cn) and a 27r-exchange component (with parameters Ci, C3 and C4). In order to treat the 
three-body correlations in isospin-asymmetric inhomogeneous nuclear many-body systems we 
assume (as done in ref.[T3]) that the relevant product of density- matrices can be represented in 
momentum space in a factorized form by T{pi, qi) r(p2, ^2) r(P3, —<li — ^2)- Such a factorization 
ansatz respects by construction the correct nuclear matter limit, but it involves approximations 
in comparison to more sophisticated treatments outlined in section 4 of ref . [TB] • Actually, the 
present approach is similar to the method DME-I introduced in ref. [18]. In comparison to 
ref.p^ the diagrammatic calculation of the isovector terms gets essentially modified only by 
relative isospin factors occurring at various places. However, their pattern is rather complex 
and therefore it is preferable to write out each (non-vanishing) contribution individually. We 
give for each diagram only the final result omitting all technical details related to extensive 
algebraic manipulations, expansions, and solving elementary integrals. 

4.1 C£;-term 

The three-body contribution from the contact interaction is represented by the left diagram in 
Fig. 1. One finds a contribution to the asymmetry energy: 

which depends quadratically on the density p = 2kj/3n'^ and is equal with opposite sign to the 
contribution to the energy per particle E{p). This property follows from the form PpPn{Pp + Pn) 
of the underlying energy density as it is determined by the Pauli exclusion principle and the 
symmetry under p -H- n exchange. Due to the momentum-independence of the three-body 
contact interaction the contributions to the other strength functions Gr,d,so,j{p) vanish. 

4.2 CD-term 

Next, we consider the three-body contributions from the Ivr-exchange component of the chiral 
3N- interaction as represented by the right diagram in Fig. 1. Putting in three medium insertions 






Figure 2: Three-body Hartree and Fock diagrams related to the chiral 27r-exchange three- nucleon 
interaction. 



one finds the following analytical expressions: 
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with the abbreviation u = kf/iriT^. Note that there is no contribution to the isovector spin-orbit 
coupling strength Gso{p), essentially because the Ivr-exchange does not generate any. 



4.3 Hartree diagram proportional to ci 3 

We continue with the three-body contributions from the 27r-exchange Hartree diagram shown 
in the left part of Fig. 2. Again putting in three medium insertions one derives the following 
analytical results: 
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which depend only on the two isoscalar coupling constants Ci and C3. The isovectorial (spin- 
dependent) C4- vertex gets eliminated by a vanishing spin-trace (over the left nucleon ring). 
The vanishing contributions to Gd{p) and Gso{p) from the 27r-exchange three-body Hartree 
diagram are particularly remarkable, in view of the fact that their isoscalar counterparts (-fd(p) 
and Fso{p) in eqs. (24,25) of ref. [13]) are quite sizeable. The actual calculation shows that the 
isospin-structure of the ci,3-vertex excludes the desired coupling of the gradient V/cp — Vfc„ to 
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4.4 Fock diagram proportional to ci 3 4 

Finally, there are the three-body contributions from the 27r-exchange Fock diagram shown in 
the right part of Fig. 2. For this diagram the occurring integrals over three Fermi spheres cannot 
be solved analytically in all cases. After a somewhat tedious calculation of the separate pieces 
proportional to Ci, C3 and C4 one finds the following results for the Fock contributions to the 
strength functions: 
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A double-index notation has been introduced for partial derivatives multiplied by powers of the 
variables x and u: 

which applies in the same way to the functions Gsij{x,u) and GTij{x,u). 
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+ — (1 + uy{A^u^ - 3) - 6(17m^ + 13m^ + 7u^ + 11) 
x^ 

+ ^(l + My-|(l + MT(3 + lln- 



+ 



7r2(8AM)4 

512m5 
4 



IOm^ + 95m2 - 79 + 16m~2 ln(l + Au^ 



15 12 

X ln2(l + 4m2) + 



2185m3 181m 
+ 



4 
119 3 

16m 2m^ 



3 3 + 16m2 



arctan 2m 
48m^ 



M-" 



16m7 
173m 9m' 



M-" 



48 



ln(l + 4m2 



ru 

/ d 
Jo 



, ,3L2 



-4(i + «Y + 4(i + nY(3n2-i) 



u 



2\2 



X^ 



x^ 



X (7 + 14m2 + 23m*) + 4x'(3m'' + 5m* - 7m' - 9) - x 



+x^(26m2 - 3m^ - 51) + x*^(2m2 - 22) 



+ 



L 



u 



15 



x^ 



(l+M 



2\5 



1 



+—(1 + m2)'(3 - 49m2) + 18(1 + 3m2)(1 + u 
x^ 

15 



2\2 



-2^^^+"')*+;j(i+" 



2n2 



(3 + 11m2)}| 



(32) 



with L{x,u) given in eq.(29). A good check of all formulas collected in this section is provided 
by their Taylor-expansion in kf. Despite the superficial opposite appearance the leading term 
in the fcj-expansion is fc^. In several cases it is even a higher power of kf. The full Taylor series 
in kf has however a small radius of convergence kf < m.,^/2, corresponding to tiny densities p < 
0.003 fm~3. Let us also take the occasion to correct the expression for Gso{p) written in eq.(27) 
of ref.y^. The correct expression is obtained by inserting into eq.(31) the parameters Ci = 0, 
C3 = —g\/2A and C4 = g\/4:A. The term omitted in ref.[TS] vanishes (accidentally) for the 
original density-matrix expansion of Negele and Vautherin [9] but not for the improved density- 
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Figure 3: Contributions to the asymmetry energy A{p) of nuclear matter. 



matrix expansion of Gebremariam, Duguet and Bogner [T^. The numerical consequences of 
this correction are insignificant since (at half nuclear matter density po/2 = 0.08 fm~^) the 
already small contribution Gsoipof^) = 6.19MeVfm^ gets just further reduced to Gsoipof^) = 
2.63MeVfm^ 



5 Results and discussion 

In this section we present and discuss our numerical results obtained by summing the series of 
two- and three-body contributions given in sections 3 and 4. The physical input parameters are: 
gA = 1-3 (nucleon axial vector coupling constant), /,r = 92.4 MeV (pion decay constant) and 
ttItt = 138 MeV (average pion mass). We use consistently the same parameters pertinent to the 
chiral three-nucleon interaction: ce = —0.625, cd = —2.06, A^ = 700 MeV, Ci = — 0.76GeV~^, 
C3 = —4.78 GeV~^ and C4 = 3.96 GeV~^, as in our previous work [T3] on the isoscalar part of the 
nuclear energy density functional. Let us remind that the low-energy constants ce = —0.625 
and Cd = —2.06 have been obtained in refs.[15[ [16] by fitting them simultaneously (after fixing 
A^ = 700 MeV) to the binding energies of ^H and ^He using the low-momentum NN-interaction 
l^ow-k [7] at a cutoff scale of A = 414 MeV. 

Fig. 3 shows the contributions to the asymmetry energy A{p) of infinite spin-saturated nu- 
clear matter for densities up to p = 0.2 fm~^. The dash-dotted and dashed line give the two-body 
and three-body contributions to this quantity. Their sum, the total interaction contribution, 
is shown by the full line in Fig. 3. In the Hartree-Fock approximation the asymmetry en- 
ergy A{p) is completed by adding the (relativistically improved) kinetic energy contribution 
-4kin(p) = kj/6M - kj/12M^, with M = 939 MeV the (free) nucleon mass. Adding these 
three pieces together, one obtains for the asymmetry energy at nuclear matter saturation den- 
sity po = 0.16 fm^^ the value A{po) = 26.5 MeV. This is compatible with the empirical values 
A{pq) = (35 ± 2) MeV extracted in extensive fits of nuclide masses in refs. [20t [2T]. For compar- 
ison, a recent microscopic estimate of the asymmetry energy in a relativistic mean-field model 
(constrained by some specific properties of certain nuclei) gave the value A{pq) = (34 ± 2) MeV 
[22] ■ Note that about 1/3 of the empirical value A{pq) is provided by the kinetic energy: 
^kin(Po) = 11.8 MeV. Another quantity of interest is the slope of the asymmetry energy at 
saturation density. We find for the slope parameter L = 3poA'{po) = 76 MeV which is again 
compatible with the value L ~ 100 MeV quoted in ref.[2T]. It is remarkable that the Hartree- 
Fock approximation works already reasonably well for the asymmetry energy A{p) when giving 
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Figure 4: Contributions to the strength function Gr{p) versus the nuclear density p 



results that are about 20% smaller than empirical determinations. In contrast to this, the 
Hartree-Fock approximation gives a much too shallow binding minimum in the equation of 
state E{p) of isospin-symmetric nuclear matter (see Fig. 7 in ref. [13]), and second order cor- 
rections are very important in order to converge eventually to the empirical saturation point 
[T^ [23]. The role of three- nucleon forces is also different for both quantities. On the one hand 
side repulsive three-body effects are essential in order to achieve saturation of nuclear matter 
but they do contribute little to the asymmetry energy A{po), in the present calculation merely 
4.0 MeV. It is also interesting to remind that the one-pion exchange alone produces a negative 
contribution to the asymmetry energy: 



Mp) 



(Itt) 



g\ml 



(47r/. 



^)2lV3 8m/ ^ ' 3 2J ' 



(33) 



which amounts to —4.4 MeV at normal nuclear matter density po = 0.16 fm~^. 

Fig. 4 shows the contributions to the strength function Gr{p)- One observes that the (neg- 
ative) two-body contributions get somewhat reduced in size by the (positive) three-body cor- 
rections. For orientation we have included in Fig. 4 also the one-pion exchange contribution to 
Gr{p) as given by the expression: 



Grip) 



(iTT) 



3(47r^) 



^!L_^_^ln(l + 4M2) 
1 + 4m2 4u ^ 



(34) 



with u = kf/rriT^. One recognizes that this relatively small contribution becomes almost inde- 
pendent of density for p > 0.02 fm~'^. Note that according its construction the Gr{p) term in 
the nuclear energy density functional eq.(2) splits the effective (in- medium) masses of protons 
and neutrons in linear proportion to a (local) isospin-asymmetry {pp — Pn)/p- By comparison 
with the isoscalar strength function Ft(p) (see Fig. 8 in ref. [13]) one concludes that the isovector 
strength function Gr{p) is suppressed by about a factor 5 and of opposite sign. 

Next, we show in Fig. 5 the strength function Gv(p) of the isovector surface term (Vpp — 
VpnY- The three-body contribution to this quantity is negligible. Furthermore, one ob- 
serves from Fig. 5 that both components Gd{p) and Gr{p)/4:p (see eq.(3)) are of equal im- 
portance for the strength function G^{p). The pronounced decrease at very low densities 
p < 0.01 fm~'^ is caused by the l7r-exchange and has also been observed in other calcula- 
tions ^ilj. In phenomenological Skyrme parameterizations the strength function G^{p)'^^^^ — 
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Figure 5: Strength function G^{p) of the isovector surface term (Vpp — Vp„) versus the nuclear 
density p. 
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Figure 6: Strength function Gso{p) of the isovector spin-orbit couphng term iy Pp—V Pn)-{Jp—Jn) 
versus the nuclear density p. 

— [3ti(2a;i + 1) + ^2(2x2 + l)]/64 is a constant, whose value is however not well determined. 
Taking the modern Sly forces [3] as a guideline one obtains the band G\j = — (11 ± 5) MeVfm^ 
which covers well the results of present (microscopic) calculation. For comparison the strength 
function F\/{p) of the isoscalar surface term (Vp)^ is about one order of magnitude larger and 
empirically much better determined: Fy ~ 75 MeVfm^ [3] . 

Next, we show in Fig. 6 the strength function Gso{p) of the isovector spin-orbit coupling term 
{y Pp—V pn)-{Jp—Jn)- One sees that the weak decrease of the two-body contribution with density 
p gets compensated by a small (positive) three-body contribution. The resulting total spin-orbit 
coupling strength Gso{p) comes out close to the constant value G^^^^ = Wq/A ~ SOMeVfm^ of 
Skyrme parameterizations [21 II]. At this point it should be emphasized that the isovector spin- 
orbit coupling strength in nuclei is presently not well determined. For example, no definite 
choice could be made in ref. [21] between different density-dependences (~ pp + 7p„, 7 = 0, 1, 2) 
of the neutron spin-orbit potential. 

Finally, we show in Fig. 7 the strength function Gj{p) accompanying the squared isovector 
spin-orbit density ( Jp — J„)^ in the nuclear energy density functional. The two-body and three- 
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Figure 7: Strength function Gj{p) multiplying the squared isovector spin-orbit density (Jp — Jn)'^ 
versus the nuclear density p. 

body contributions come with equal sign but exhibit an opposite density-dependence. For 
orientation we have included in Fig. 7 also the one-pion exchange contribution to Gj{p) as given 
by the expression: 



gAp) 



(1^) 



3^i 



(32m^^)%e 



4u^- 



Su 



- ln(l + Au^ 



(35) 



Apparently, the strength function Gj{p) is dominated by this unique long-range contribution 
which is also responsible for the strong density-dependence fTI] below p < 0.05 fm~^. At this 
point it should be kept in mind that the {Jp — JnY term in the nuclear energy density func- 
tional eq.(2) represents non-local Fock contributions from tensor forces etc. An outstanding 
iTT-exchange contribution to Gj{p) is therefore not surprising. 

In summary we have calculated the strength functions of isovector terms in the nuclear 
energy density functional from chiral two-and three-nucleon interactions. The results for asym- 
metry energy A{p) suggest that the Hartree-Fock approximation could work better for isovector 
quantities. Clearly, it remains a challenge to confirm this in a consistent second-order calculation 
of the complete nuclear energy density functional. 
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